/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _GENERALEQUATION_H_
#define _GENERALEQUATION_H_

#include <string>
#include <map>
#include "Tuna.hpp"
#include "Storage/SparseMatrix.hpp"
#include "Storage/Diagonal.hpp"
#include "Utils/num_utils.hpp"

namespace Tuna {

/*! 
 ***************************************************************************
 *  General class for equations.
 *  This class is a model for the concept "General Equation" defined
 *  in the analysis and design of my doctoral thesis. It is the base
 *  class to derive all discretized governing equations, using 
 *  Finite Volume Method.
 *  
 ***************************************************************************
 * \author  Luis M. de la Cruz [ Sat May  2 18:11:12 CDT 2009 ]
 ***************************************************************************
 */  
 
template<typename Teq>
class GeneralEquation {
public:
  
/*! Precision used for the numerical operations. */
  typedef typename Typeinfo<Teq>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Teq>::Dim >::tiny floatTinyArray_t; 
  typedef typename TunaArray<prec_t, Typeinfo<Teq>::Dim >::huge ScalarField; 
  typedef typename TunaArray<prec_t,1>::huge ScalarField1D;
  typedef map<BC_t, prec_t> BC_mapping;  
  typedef SparseMatrix< Diagonal< prec_t, Typeinfo<Teq>::Dim > > DiagMat;


  GeneralEquation();
  GeneralEquation(ScalarField&, DiagMat&, ScalarField&);
  GeneralEquation(ScalarField&, DiagMat&, ScalarField&, const floatTinyArray_t&);
  ~GeneralEquation() { }; 
/*!
 *  The Curiously Recursive Template Pattern (CRTP) is used.
 *  asDerived() gets an returns a reference to the derived class.
 */
  inline Teq& asDerived() { return static_cast<Teq&>(*this); }
   
//  
// ----- Setup parameters of the equation
//

  inline void setLinearSystem(DiagMat&, ScalarField &);

  inline void setDeltas(const floatTinyArray_t& deltas) {
    dx = deltas(0); 
    if( Dim > 1) dy = deltas(1); 
    if( Dim > 2) dz = deltas(2);
  }  

  template<typename Tmesh>
  inline void setNonUniformMesh(Tmesh& mesh) {
    dxyz[0].resize( mesh.getExtentDeltasVols(X) );
    dxyz[1].resize( mesh.getExtentDeltasVols(Y) );
    dxyz[2].resize( mesh.getExtentDeltasVols(Z) );
    dxyz[0] = mesh.getDeltasVols(X);
    dxyz[1] = mesh.getDeltasVols(Y);
    dxyz[2] = mesh.getDeltasVols(Z);    

    dface[0].resize( mesh.getExtentDeltasFace(X) );
    dface[1].resize( mesh.getExtentDeltasFace(Y) );
    dface[2].resize( mesh.getExtentDeltasFace(Z) );
    dface[0] = mesh.getDeltasFace(X);
    dface[1] = mesh.getDeltasFace(Y);
    dface[2] = mesh.getDeltasFace(Z);

    dx = dxyz[0](1);
    if (Dim >= 2) dy = dxyz[1](1);
    if (Dim == 3) dz = dxyz[2](1);

    cout << "\n dx = " << dx;
    cout << "\n dy = " << dy;
    cout << "\n dz = " << dz;
  }

  inline void setDeltaTime(prec_t ddt) { dt = ddt; }  
//
// ----- Setup and apply boundary conditions
//
  inline void setDirichlet(BC_t wall, prec_t w_v = 0.0);
  inline void setNeumann(BC_t wall, prec_t w_v = 0.0);  
  inline void clearDirichlet() { dirichlet.clear();}
  inline void clearNeumann() { neumann.clear();}
  inline void applyBoundaryConditions();
  inline void applyBoundaryConditions1D();
  inline void applyBoundaryConditions2D();
  inline void applyBoundaryConditions3D();
  inline void applyDirichlet1D();
  inline void applyDirichlet2D();
  inline void applyDirichlet3D();
  inline void applyNeumann1D();
  inline void applyNeumann2D();
  inline void applyNeumann3D();

  inline void applyBoundaryConditionsStagger1D();
  inline void applyBoundaryConditionsStagger2D(Axis_t);
  inline void applyBoundaryConditionsStagger3D(Axis_t);
  inline void applyDirichletStagger1D();
  inline void applyDirichletStagger2D(Axis_t);
  inline void applyDirichletStagger3D(Axis_t);

  inline void set_bi(int bbi) {bi = bbi;}
  inline void set_bj(int bbj) {bj = bbj;}
  inline void set_bk(int bbk) {bk = bbk;}
  inline void set_ei(int eei) {ei = eei;}
  inline void set_ej(int eej) {ej = eej;}
  inline void set_ek(int eek) {ek = eek;}

//
// ----- Get
//
  inline const ScalarField& get_phi() const { return phi; }
  inline ScalarField& get_phi(Range &r) { return phi(r); }
  inline const ScalarField& get_aP(Range &r) const { return aP(r); }
  inline const ScalarField& get_aE(Range &r) const { return aE(r); }
  inline const ScalarField& get_aW(Range &r) const { return aW(r); }
  inline const ScalarField& get_aN(Range &r) const { return aN(r); }
  inline const ScalarField& get_aS(Range &r) const { return aS(r); }
  inline const ScalarField& get_aF(Range &r) const { return aF(r); }
  inline const ScalarField& get_aB(Range &r) const { return aB(r); }
  inline const ScalarField& get_sp(Range &r) const { return sp(r); }
  inline int get_bi() const { return bi; }
  inline int get_ei() const { return ei; }
  inline int get_bj() const { return bj; }
  inline int get_ej() const { return ej; }
  inline int get_bk() const { return bk; }
  inline int get_ek() const { return ek; }
  inline prec_t getDx() { return dx; }
  inline prec_t getDy() { return dy; }
  inline prec_t getDz() { return dz; }
  inline prec_t getDt() { return dt; }

//
// ----- Useful functions
//
  inline double calcResidual();
  inline prec_t calcErrorL2() { return NumUtils::calcErrorL2(phi_0, phi); }
  inline prec_t calcErrorL1() { return NumUtils::calcErrorL1(phi_0, phi); }
  inline prec_t calcErrorLmax() { return NumUtils::calcErrorLmax(phi_0, phi); }

  double getResidual() { return residual; }
  void update();
  void updatePhi(ScalarField &s) { phi = s; } 
  inline bool applyBounds(int bbi, int eei)
  {
    bi = bbi; ei = eei;
  }
  inline bool applyBounds(int bbi, int eei,
			  int bbj, int eej)
  {
    bi = bbi; ei = eei; 
    bj = bbj; ej = eej; 
  } 
  inline bool applyBounds(int bbi, int eei,
			  int bbj, int eej,
			  int bbk, int eek)
  {
    bi = bbi; ei = eei; 
    bj = bbj; ej = eej; 
    bk = bbk; ek = eek; 
  }  
  void print(int rank = -1);

protected:
  int Dim;
  string name;
  prec_t dx, dy, dz, dt, residual;
  ScalarField1D dxyz[3], dface[3];  
  int bi, ei, bj, ej, bk, ek;
  BC_mapping dirichlet, neumann;
        
public:
  ScalarField phi;
  ScalarField aE, aW, aN, aS, aF, aB, aP, sp;

private:
  ScalarField phi_0;
};
  
/*!
 *  Default constructor.
 */
template<typename Teq>
GeneralEquation<Teq>::GeneralEquation ()
{
  Dim = Typeinfo<Teq>::Dim;
  dx = 1.0; dy = 1.0; dz = 1.0; dt = 1.0;
  bi = 0; ei = 1; bj = 0; ej = 1; bk = 0; ek = 1;
} 
  
/*!
 *  In this constructor \c phi is a copy of a global field variable.
 *  \c phi is used to storage the solution in each iteration of the solver.
 */
template<typename Teq>
GeneralEquation<Teq>::GeneralEquation(ScalarField &phi_global, 
				      DiagMat &matrix, ScalarField &b)
{
  Dim = Typeinfo<Teq>::Dim;
  dx = 1.0; dy = 1.0; dz = 1.0; dt = 1.0;

  // phi is a copy of a global field variable and 
  // store the new solution.
  phi.resize(phi_global.shape());
  phi = phi_global;

  // phi_0 is a reference of a global field variable
  // and stores the previous solution
  phi_0.reference(phi_global);

  bi = bj = bk = 0; ei = ej = ek = 1;

  if (Dim >= 1) {
    bi = phi.lbound(firstDim);  ei = phi.ubound(firstDim);
  } 
  if (Dim >= 2) {
    bj = phi.lbound(secondDim); ej = phi.ubound(secondDim);
  }
  if (Dim == 3) {
    bk = phi.lbound(thirdDim);  ek = phi.ubound(thirdDim);
  }
  setLinearSystem(matrix, b);
}
  
/*!
 *  In this constructor \c phi is a copy of a global field variable.
 *  \c phi is used to storage the solution in each iteration of the solver.
 */
template<typename Teq>
GeneralEquation<Teq>::GeneralEquation(ScalarField &phi_global, 
				      DiagMat &matrix, ScalarField &b,
				      const floatTinyArray_t& deltas) 
{
  Dim = Typeinfo<Teq>::Dim;
  dx = 1.0; dy = 1.0; dz = 1.0; dt = 1.0;

  // phi is a copy of a global field variable and 
  // store the new solution.
  phi.resize(phi_global.shape());
  phi = phi_global;

  // phi_0 is a reference of a global field variable
  // and stores the previous solution
  phi_0.reference(phi_global);

  bi = bj = bk = 0; ei = ej = ek = 1;

  if (Dim >= 1) {
    bi = phi.lbound(firstDim);  ei = phi.ubound(firstDim);
  } 
  if (Dim >= 2) {
    bj = phi.lbound(secondDim); ej = phi.ubound(secondDim);
  }
  if (Dim == 3) {
    bk = phi.lbound(thirdDim);  ek = phi.ubound(thirdDim);
  }
  setLinearSystem(matrix, b);
  setDeltas(deltas);
}

/*!
 *  Setup the Linear System.
 *  The arrays \c aE, \c aW, \c aN, \c aS, \c aF, \c aB, \c aP and \c sp are 
 *  references to global matrices and vectors used to storage the linear system.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::setLinearSystem
(DiagMat &matrix, ScalarField &b )
{
  aP.reference(matrix.store(0)); // main diagonal
  aE.reference(matrix.store(1)); // upper diagonal
  aW.reference(matrix.store(2)); // lower diagonal
  if(Dim >= 2) {
    aN.reference(matrix.store(3)); // second upper diagonal
    aS.reference(matrix.store(4)); // second lower diagonal
  }
  if(Dim == 3) {
    aF.reference(matrix.store(5)); // third upper diagonal
    aB.reference(matrix.store(6)); // third lower diagonal
  }
  sp.reference(b); // right-hand side 
}

template<typename Teq>
void GeneralEquation<Teq>::print(int rank) 
{
  if (rank == -1)
    std::cout << "\n +----- TUNA info -----+"
	      << "\n | GeneralEquation< " << name << "  > "
	      << "\n +-----+";
  else
    std::cout << "\n +----- TUNA info -----+"
	      << "\n | Proc (" << rank << ") : GeneralEquation< " << name << "  > "
	      << "\n +-----+";

  std::cout << "\n | bi = " << bi << "\t ei = " << ei
	    << "\n | bj = " << bj << "\t ej = " << ej
	    << "\n | bk = " << bk << "\t ek = " << ek
	    << "\n | dx = " << dx << "\t dy = " << dy
	    << "\n | dz = " << dz << "\t dt = " << dt    
	    << "\n +-----+"
	    << "\n | Boundary conditions "
	    << "\n +-----+";

    std::string bc_name[6] = {"TOP", "BOT", "LEFT", "RIGHT", "FRONT", "BACK"};
    
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
      std::cout<< "\n | Dirichlet: Wall = " << bc_name[pos->first];
    }
    for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
      std::cout<< "\n | Neumann: Wall = " << bc_name[pos->first];
    }
    asDerived().printInfo();
}

/*!
 *  Setup Dirichlet boundary conditions on walls.
 *  This function just moves the beginnings (bi,bj,bk) and endings (ei,ej,ek)
 *  of the Scalarfield arrays, and store the name of the wall where the 
 *  condition was applied in the \c dirichlet container (stl::map).
 */
template<typename Teq>
inline void GeneralEquation<Teq>::setDirichlet(BC_t wall, prec_t w_v)
{
    dirichlet[wall] = w_v;
    switch(wall) {
    case RIGHT_WALL : ei--; break;
    case LEFT_WALL  : bi++; break;
    case TOP_WALL   : ej--; break;
    case BOTTOM_WALL: bj++; break;
    case FRONT_WALL : ek--; break;
    case BACK_WALL  : bk++; break;
    default:
	cout << " GeneralEquation: setDirichlet: wall " 
	     << wall << " not defined "<< endl;
    } 
}

/*!
 *  Setup Neumann boundary conditions on walls.
 *  This function just moves the beginnings (bi,bj,bk) and endings (ei,ej,ek)
 *  of the Scalarfield arrays, and store the name of the wall where the 
 *  condition was applied in the \c neumann container (stl::map).
 */
template<typename Teq>
inline void GeneralEquation<Teq>::setNeumann(BC_t wall, prec_t w_v)
{
    neumann[wall] = w_v;
    switch(wall) {
    case RIGHT_WALL : ei--; break;
    case LEFT_WALL  : bi++; break;
    case TOP_WALL   : ej--; break;
    case BOTTOM_WALL: bj++; break;
    case FRONT_WALL : ek--; break;
    case BACK_WALL  : bk++; break;
    default:
	cout << " GeneralEquation: setNeumann: wall " 
	     << wall << " not defined "<< endl;
    }
}

template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditions() {
  if(Dim == 1) {
    applyDirichlet1D();
    applyNeumann1D();
  } else if (Dim == 2) {
    applyDirichlet2D();
    applyNeumann2D();
  } else if (Dim == 3) {  
    applyDirichlet3D();
    applyNeumann3D(); 
  }
}
/*!
 *  Application of Dirichlet and Neumann boundary conditions in 1D. 
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditions1D() {
    applyDirichlet1D();
    applyNeumann1D();
}

/*!
 *  Application of Dirichlet and Neumann boundary conditions in 2D. 
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditions2D() {
    applyDirichlet2D();
    applyNeumann2D();
}

/*!
 *  Application of Dirichlet and Neumann boundary conditions in 3D. 
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditions3D() {
    applyDirichlet3D();
    applyNeumann3D();
}
  

/************************************************************
 *  The same as above but for staggered meshes
/************************************************************/
template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditionsStagger1D() {
    applyDirichletStagger1D();
    applyNeumann1D();
}

template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditionsStagger2D(Axis_t axis) {
    applyDirichletStagger2D(axis);
    applyNeumann2D();
}

template<typename Teq>
inline void GeneralEquation<Teq>::applyBoundaryConditionsStagger3D(Axis_t axis) {
    applyDirichletStagger3D(axis);
    applyNeumann3D();
}


/*!
 *  Dirichlet boundary condition in 1D. 
 *  This function applies the Dirichlet boundary conditions iterating on the 
 *  \c dirichlet container that was defined by the setDirichlet() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichlet1D()
{
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei) += aE(ei);
	    sp(ei) += 2 * aE(ei) * phi(ei+1);
	    aE(ei) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi) += aW(bi);
	    sp(bi) += 2 * aW(bi) * phi(bi-1);
	    aW(bi) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "<< endl;
	}
    }    
}

/************************************************************
 *  The same as above but for staggered meshes
/************************************************************/
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichletStagger1D()
{
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    sp(ei) += aE(ei) * phi(ei+1);
	    aE(ei) = 0.0;
	    break;
	case LEFT_WALL:
	    sp(bi) += aW(bi) * phi(bi-1);
	    aW(bi) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "<< endl;
	}
    }    
}

/*!
 *  Dirichlet boundary condition in 2D. 
 *  This function applies the Dirichlet boundary conditions iterating on the 
 *  \c dirichlet container that was defined by the setDirichlet() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichlet2D() {
    Range I(bi,ei), J(bj, ej);
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei, J) += aE(ei, J);
	    sp(ei, J) += 2 * aE(ei, J) * phi(ei+1, J);
	    aE(ei, J) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi, J) += aW(bi, J);		
	    sp(bi, J) += 2 * aW(bi, J) * phi(bi-1, J);
	    aW(bi, J) = 0.0;
	    break;
	case TOP_WALL:
	    aP(I, ej) += aN(I, ej);
	    sp(I, ej) += 2 * aN(I, ej) * phi(I, ej+1);
	    aN(I, ej) = 0.0;
	    break;
	case BOTTOM_WALL:
	    aP(I, bj) += aS(I, bj);
	    sp(I, bj) += 2 * aS(I, bj) * phi(I, bj-1);
	    aS(I, bj) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "<< endl;
	}
    }

}

/************************************************************
 *  The same as above but for staggered meshes
/************************************************************/
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichletStagger2D(Axis_t axis) {
    Range I(bi,ei), J(bj, ej);
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	  if(axis == 0) {
	    sp(ei, J) += aE(ei, J) * phi(ei+1, J);
	    aE(ei, J) = 0.0;
	  }
	  else {
	    aP(ei, J) += aE(ei, J);
	    sp(ei, J) += 2 * aE(ei, J) * phi(ei+1, J);
	    aE(ei, J) = 0.0;
	  }
	  break;
	case LEFT_WALL:
	  if(axis == 0) {
	    sp(bi, J) += aW(bi, J) * phi(bi-1, J);
	    aW(bi, J) = 0.0;
	  } 
	  else {
	    aP(bi, J) += aW(bi, J);		
	    sp(bi, J) += 2 * aW(bi, J) * phi(bi-1, J);
	    aW(bi, J) = 0.0;	    
	  }
	  break;
	case TOP_WALL:
	  if(axis == 1) {
	    sp(I, ej) += aN(I, ej) * phi(I, ej+1);
	    aN(I, ej) = 0.0;
	  } 
	  else {
	    aP(I, ej) += aN(I, ej);
	    sp(I, ej) += 2 * aN(I, ej) * phi(I, ej+1);
	    aN(I, ej) = 0.0;	    
	  }
	  break;
	case BOTTOM_WALL:
	  if(axis == 1) {
	    sp(I, bj) += aS(I, bj) * phi(I, bj-1);
	    aS(I, bj) = 0.0;
	  }
	  else {
	    aP(I, bj) += aS(I, bj);
	    sp(I, bj) += 2 * aS(I, bj) * phi(I, bj-1);
	    aS(I, bj) = 0.0;
	  }
	  break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "<< endl;
	}
    }
}

/*!
 *  Dirichlet boundary condition in 3D. 
 *  This function applies the Dirichlet boundary conditions iterating on the 
 *  \c dirichlet container that was defined by the setDirichlet() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichlet3D()
{
    Range I(bi,ei), J(bj, ej), K(bk, ek);
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei,J,K) += aE(ei, J, K);
	    sp(ei,J,K) += 2 * aE(ei,J,K) * phi(ei+1,J,K);
	    aE(ei,J,K) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi,J,K) += aW(bi, J, K);
	    sp(bi,J,K) += 2 * aW(bi,J,K) * phi(bi-1,J,K);
	    aW(bi,J,K) = 0.0;
	    break;
	case TOP_WALL:
	    aP(I,ej,K) += aN(I, ej, K);
	    sp(I,ej,K) += 2 * aN(I,ej,K) * phi(I,ej+1,K);
	    aN(I,ej,K) = 0.0;
	    break;
	case BOTTOM_WALL:
	    aP(I,bj,K) += aS(I, bj, K);
	    sp(I,bj,K) += 2 * aS(I,bj,K) * phi(I,bj-1,K);
	    aS(I,bj,K) = 0.0;
	    break;
	case FRONT_WALL:
	    aP(I,J,ek) += aF(I, J, ek);
	    sp(I,J,ek) += 2 * aF(I,J,ek) * phi(I,J,ek+1); 
	    aF(I,J,ek) = 0.0;
	    break;
	case BACK_WALL:
	    aP(I,J,bk) += aB(I, J, bk);
	    sp(I,J,bk) += 2 * aB(I,J,bk) * phi(I,J,bk-1);
	    aB(I,J,bk) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "
		 << name << endl;
	}
    }
}

/************************************************************
 *  The same as above but for staggered meshes
/************************************************************/
template<typename Teq>
inline void GeneralEquation<Teq>::applyDirichletStagger3D(Axis_t axis)
{
    Range I(bi,ei), J(bj, ej), K(bk, ek);
    typename BC_mapping::iterator pos;
    for(pos = dirichlet.begin(); pos !=dirichlet.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	  if(axis == 0) {
	    sp(ei,J,K) += aE(ei,J,K) * phi(ei+1,J,K);
	    aE(ei,J,K) = 0.0;
	  } 
	  else {
	    aP(ei,J,K) += aE(ei, J, K);
	    sp(ei,J,K) += 2 * aE(ei,J,K) * phi(ei+1,J,K);
	    aE(ei,J,K) = 0.0;
	  }
	  break;
	case LEFT_WALL:
	  if(axis == 0) {
	    sp(bi,J,K) += aW(bi,J,K) * phi(bi-1,J,K);
	    aW(bi,J,K) = 0.0;	    
	  }
	  else {
	    aP(bi,J,K) += aW(bi, J, K);
	    sp(bi,J,K) += 2 * aW(bi,J,K) * phi(bi-1,J,K);
	    aW(bi,J,K) = 0.0;
	  }
	  break;
	case TOP_WALL:
	  if(axis == 1) {
	    sp(I,ej,K) += aN(I,ej,K) * phi(I,ej+1,K);
	    aN(I,ej,K) = 0.0;
	  }
	  else {
	    aP(I,ej,K) += aN(I, ej, K);
	    sp(I,ej,K) += 2 * aN(I,ej,K) * phi(I,ej+1,K);
	    aN(I,ej,K) = 0.0;
	  }
	  break;
	case BOTTOM_WALL:
	  if(axis == 1) {
	    sp(I,bj,K) += aS(I,bj,K) * phi(I,bj-1,K);
	    aS(I,bj,K) = 0.0;
	  }
	  else {
	    aP(I,bj,K) += aS(I, bj, K);
	    sp(I,bj,K) += 2 * aS(I,bj,K) * phi(I,bj-1,K);
	    aS(I,bj,K) = 0.0;
	  }
	  break;
	case FRONT_WALL:
	  if(axis == 2) {
	    sp(I,J,ek) += aF(I,J,ek) * phi(I,J,ek+1); 
	    aF(I,J,ek) = 0.0;	    
	  }
	  else {
	    aP(I,J,ek) += aF(I, J, ek);
	    sp(I,J,ek) += 2 * aF(I,J,ek) * phi(I,J,ek+1); 
	    aF(I,J,ek) = 0.0;
	  }
	  break;
	case BACK_WALL:
	  if(axis == 2) {
	    sp(I,J,bk) += aB(I,J,bk) * phi(I,J,bk-1);
	    aB(I,J,bk) = 0.0;	    
	  } 
	  else {
	    aP(I,J,bk) += aB(I, J, bk);
	    sp(I,J,bk) += 2 * aB(I,J,bk) * phi(I,J,bk-1);
	    aB(I,J,bk) = 0.0;
	  }
	  break;
	default:
	    cout << " GeneralEquation: applyDirichlet: wall " 
		 << pos->first << " not defined "
		 << name << endl;
	}
    }
}

/*!
 *  Neumann boundary condition in 1D. 
 *  This function applies the Neumann boundary conditions iterating on the 
 *  \c neumann container that was defined by the setNeumann() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyNeumann1D()
{ 
    typename BC_mapping::iterator pos;
    for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei) -= aE(ei);
	    sp(ei) += aE(ei) * dx * neumann[RIGHT_WALL];
	    aE(ei) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi) -= aW(bi);
	    sp(bi) -= aW(bi) * dx * neumann[LEFT_WALL];
	    aW(bi) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyNeumann : wall " 
		 << pos->first << " not defined "<< endl;
	}
    }
}


/*!
 *  Neumann boundary condition in 2D. 
 *  This function applies the Neumann boundary conditions iterating on the 
 *  \c neumann container that was defined by the setNeumann() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyNeumann2D()
{
  Range I(bi,ei), J(bj, ej);
  //    Range I(bi-1,ei+1), J(bj-1, ej+1);
    typename BC_mapping::iterator pos;
    for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei, J) -= aE(ei, J);
	    sp(ei, J) += aE(ei, J) * dx * neumann[RIGHT_WALL];
	    aE(ei, J) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi, J) -= aW(bi, J);
	    sp(bi, J) -= aW(bi, J) * dx * neumann[LEFT_WALL];
	    aW(bi, J) = 0.0;
	    break;
	case TOP_WALL:
	    aP(I, ej) -= aN(I, ej);
	    sp(I, ej) += aN(I, ej) * dy * neumann[TOP_WALL];
	    aN(I, ej) = 0.0;
	    break;
	case BOTTOM_WALL:
	    aP(I, bj) -= aS(I, bj);
	    sp(I, bj) -= aS(I, bj) * dy * neumann[BOTTOM_WALL];
	    aS(I, bj) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyNeumann : wall " 
		 << pos->first << " not defined "<< endl;
	}
    }
}

/*!
 *  Neumann boundary condition in 3D. 
 *  This function applies the Neumann boundary conditions iterating on the 
 *  \c neumann container that was defined by the setNeumann() function. The
 *  coeficients aE, aW, aN, aS, aN, aF, aB, aP and sp (source), are modified 
 *  according to these conditions.
 */
template<typename Teq>
inline void GeneralEquation<Teq>::applyNeumann3D() {
    Range I(bi,ei), J(bj, ej), K(bk, ek);
    typename BC_mapping::iterator pos;
    for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	    aP(ei, J, K) -= aE(ei, J, K);
	    sp(ei, J, K) += aE(ei, J, K) * dx * neumann[RIGHT_WALL];
	    aE(ei, J, K) = 0.0;
	    break;
	case LEFT_WALL:
	    aP(bi, J, K) -= aW(bi, J, K);
	    sp(bi, J, K) -= aW(bi, J, K) * dx * neumann[LEFT_WALL];
	    aW(bi, J, K) = 0.0;
	    break;
	case TOP_WALL:
	    aP(I, ej, K) -= aN(I, ej, K);
	    sp(I, ej, K) += aN(I, ej, K) * dy * neumann[TOP_WALL];
	    aN(I, ej, K) = 0.0;
	    break;
	case BOTTOM_WALL:
	    aP(I, bj, K) -= aS(I, bj, K);
	    sp(I, bj, K) -= aS(I, bj, K) * dy * neumann[BOTTOM_WALL];
	    aS(I, bj, K) = 0.0;
	    break;
	case FRONT_WALL:
	    aP(I, J, ek) -= aF(I, J, ek);
	    sp(I, J, ek) += aF(I, J, ek) * dz * neumann[FRONT_WALL];
	    aF(I, J, ek) = 0.0;
	    break;
	case BACK_WALL:
	    aP(I, J, bk) -= aB(I, J, bk);
	    sp(I, J, bk) -= aB(I, J, bk) * dz * neumann[BACK_WALL];
	    aB(I, J, bk) = 0.0;
	    break;
	default:
	    cout << " GeneralEquation: applyNeumann : wall " 
		 << pos->first << " not defined "<< endl;
	}
    }
}

template<typename Teq>
inline void GeneralEquation<Teq>::update()
{

  // First we apply the Neumann conditions, if any...
    typename BC_mapping::iterator pos;
    Range I(bi-1,ei+1), J(bj-1, ej+1), K(bk-1, ek+1);

    if (Dim == 1) {
      for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	  phi(ei+1) = phi(ei) + neumann[RIGHT_WALL] * dx * 0.5;
	  break;
	case LEFT_WALL:
	  phi(bi-1) = phi(bi) - neumann[LEFT_WALL] * dx * 0.5;
	  break;
	default:
	  cout << " GeneralEquation: applyNeumann : wall " 
	       << pos->first << " not defined "<< endl;
	}
      }
    }

    if (Dim == 2) {
      for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	  phi(ei+1, J) = phi(ei, J) + neumann[RIGHT_WALL] * dx * 0.5;
	  break;
	case LEFT_WALL:
	  phi(bi-1, J) = phi(bi, J) - neumann[LEFT_WALL] * dx * 0.5;
	  break;
	case TOP_WALL:
	  phi(I, ej+1) = phi(I, ej) + neumann[TOP_WALL] * dy * 0.5;
	  break;
	case BOTTOM_WALL:
	  phi(I, bj-1) = phi(I, bj) - neumann[BOTTOM_WALL] * dy * 0.5;
	  break;
	default:
	    cout << " GeneralEquation: applyNeumann : wall " 
		 << pos->first << " not defined "<< endl;
	}
      }
    }

    if (Dim == 3) {
      for(pos = neumann.begin(); pos !=neumann.end(); ++pos) {
	switch(pos->first) {
	case RIGHT_WALL:
	  phi(ei+1, J, K) = phi(ei, J, K) + neumann[RIGHT_WALL] * dx * 0.5;
	  break;
	case LEFT_WALL:
	  phi(bi-1, J, K) = phi(bi, J, K) - neumann[LEFT_WALL] * dx * 0.5;
	  break;
	case TOP_WALL:
	  phi(I, ej+1, K) = phi(I, ej, K) + neumann[TOP_WALL] * dy * 0.5;
	  break;
	case BOTTOM_WALL:
	  phi(I, bj-1, K) = phi(I, bj, K) - neumann[BOTTOM_WALL] * dy * 0.5;
	  break;
	case FRONT_WALL:
	  phi(I, J, ek+1) = phi(I, J, ek) + neumann[FRONT_WALL] * dz * 0.5;
	  break;
	case BACK_WALL:
	  phi(I, J, bk-1) = phi(I, J, bk) - neumann[BACK_WALL] * dz * 0.5;
	  break;
	default:
	    cout << " GeneralEquation: applyNeumann : wall " 
		 << pos->first << " not defined "<< endl;
	}
      }
    }

    // update phi_0 global using the solution stored in phi.   
    phi_0 = phi;
}

template<typename Teq>
double GeneralEquation<Teq>::calcResidual() 
{
    Range I(bi, ei), J(bj, ej), K(bk, ek);	
    static ScalarField acum(aP.shape());
    acum = 0.0;
    
    int num_grid_points = (ei - bi) * (ej - bj) * (ek -bk);
    double NINV = 1.0 / num_grid_points;

    if ( Dim == 1 ) {
	acum(I) = aE(I) * phi(I+1) + aW(I) * phi(I-1) +
	    sp(I) - aP(I) * phi(I);
    }
    
    if ( Dim == 2 ) {
	acum(I,J) = aE(I,J) * phi(I+1,J) + aW(I,J) * phi(I-1,J) +
	    aN(I,J) * phi(I,J+1) + aS(I,J) * phi(I,J-1) +
	    sp(I,J) - aP(I,J) * phi(I,J);
    }

    if ( Dim == 3 ) {
	acum(I,J,K) = aE(I,J,K) * phi(I+1,J,K) + aW(I,J,K) * phi(I-1,J,K)
	    + aN(I,J,K) * phi(I,J+1,K) + aS(I,J,K) * phi(I,J-1,K) 
	    + aF(I,J,K) * phi(I,J,K+1) + aB(I,J,K) * phi(I,J,K-1);
	acum(I,J,K) += sp(I,J,K) - aP(I,J,K) * phi(I,J,K);
    }    
// L2-Norm
    residual = sqrt( NINV * sum(acum * acum) );
//	residual = NINV * sqrt( sum(acum * acum) );
// L1-Norm
//	residual = NINV * sum( abs(acum) );
// Lmax-Norm
//	residual = max ( abs(acum) );
    return residual;
}


} // namespace Tuna

#endif //_GENERALEQUATION_H_











